Son muchos los paquetes empleados en estos análisis. Puedes consultar en el ChatGPT qué hace cada uno. Considera un aspecto también importante: algunas funciones escritas por mí se cargan con source_url y source; dentro de algunas de dichas funciones, también se cargan paquetes adicionales.
library(vegan)
library(sf)
library(tidyverse)
library(tmap)
library(kableExtra)
library(broom)
library(cluster)
library(gclus)
library(pvclust)
library(foreach)
library(leaps)
library(caret)
library(RColorBrewer)
library(indicspecies)
library(dendextend)
library(adespatial)
library(SpadeR)
library(iNEXT)
library(GGally)
library(vegetarian)
r <- 'R/'
gh_content <- 'https://raw.githubusercontent.com/'
gh_zonal_stats <- paste0(gh_content,
'geofis/zonal-statistics/0b2e95aaee87bf326cf132d28f4bd15220bb4ec7/out/')
repo_analisis <- 'biogeografia-master/scripts-de-analisis-BCI/master'
repo_sem202202 <- 'biogeografia-202202/material-de-apoyo/master/practicas/'
devtools::source_url(paste0(gh_content, repo_analisis, '/biodata/funciones.R'))
devtools::source_url(paste0(gh_content, repo_sem202202, 'train.R'))
devtools::source_url(paste0(gh_content, repo_sem202202, 'funciones.R'))
fuentes_practica <- 'fuentes/practica-03/'
source(paste0(r, 'funciones.R'))
umbral_alfa <- 0.05
mc <- read.csv(paste0(fuentes_practica, 'matriz-comunidad-', params$estudiante, '.csv'))[,-1]
mc %>% estilo_kable(
titulo = 'Matriz de comunidad',
nombres_filas = T, alinear = 'r')
| sp01 | sp02 | sp03 | sp04 | sp05 | sp06 | sp07 | sp08 | sp09 | sp10 | sp11 | sp12 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 0 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 0 |
| 2 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 1 | 1 |
| 3 | 0 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 |
| 4 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 0 |
| 5 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 |
| 6 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 1 |
| 7 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 1 | 0 | 1 |
| 8 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 |
| 9 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 |
| 10 | 1 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 1 | 0 | 0 |
| 11 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 1 |
| 12 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 |
| 13 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 1 |
| 14 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 1 | 1 |
| 15 | 1 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 1 |
data.frame(Especies = sort(names(mc))) %>%
estilo_kable(titulo = 'Lista de especies', cubre_anchura = F, alinear = 'c') %>%
column_spec(column = 1, width = "15em")
| Especies |
|---|
| sp01 |
| sp02 |
| sp03 |
| sp04 |
| sp05 |
| sp06 |
| sp07 |
| sp08 |
| sp09 |
| sp10 |
| sp11 |
| sp12 |
data.frame(`Número de sitios donde fue reportada la especie` = sort(colSums(mc), decreasing = T),
check.names = F) %>%
rownames_to_column('Especie') %>%
estilo_kable(
titulo = 'Número de sitios en los que está presente cada especie (orden descendente por número de sitios)',
nombres_filas = F, alinear = 'cr')
| Especie | Número de sitios donde fue reportada la especie |
|---|---|
| sp04 | 11 |
| sp06 | 11 |
| sp07 | 11 |
| sp09 | 11 |
| sp10 | 11 |
| sp03 | 10 |
| sp11 | 10 |
| sp12 | 10 |
| sp01 | 9 |
| sp05 | 9 |
| sp02 | 8 |
| sp08 | 8 |
data.frame(`Riqueza por sitios` = rowSums(mc),
check.names = F) %>% rownames_to_column('Sitio') %>%
arrange(desc(`Riqueza por sitios`)) %>%
estilo_kable(
titulo = 'Riqueza por sitios (orden descendente por riqueza)',
nombres_filas = F, alinear = 'cr')
| Sitio | Riqueza por sitios |
|---|---|
| 9 | 11 |
| 11 | 10 |
| 1 | 9 |
| 7 | 9 |
| 8 | 9 |
| 14 | 9 |
| 2 | 8 |
| 3 | 8 |
| 4 | 8 |
| 5 | 8 |
| 13 | 8 |
| 15 | 7 |
| 6 | 6 |
| 10 | 6 |
| 12 | 3 |
La matriz de comunidad analizada se compone de 15 sitios y 12 especies, donde el/los sitio/s más ricos es/son 9. La/s especie/s más común/es es/son sp04, sp06, sp07, sp09 y sp10 y la/s más rara/s es/son sp02 y sp08. El siguiente gráfico de mosaicos muestra la distribución de las especies según sitios.
grafico_mosaico <- crear_grafico_mosaico_de_mc(mc, tam_rotulo = 12) + xlab('Sitios') + ylab('Especie')
grafico_mosaico
FIGURA 1.1: Distribución de las especies según sitios
Este paso es importante, lo explico aquí
mc_t <- decostand(mc, 'hellinger') #Hellinger, funciona con datos de presencia/ausencia
mc_t %>% estilo_kable(titulo = 'Matriz de comunidad transformada',
nombres_filas = T, alinear = 'r')
| sp01 | sp02 | sp03 | sp04 | sp05 | sp06 | sp07 | sp08 | sp09 | sp10 | sp11 | sp12 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0.33 | 0.33 | 0.33 | 0.00 | 0.33 | 0.33 | 0.00 | 0.33 | 0.33 | 0.33 | 0.33 | 0.00 |
| 2 | 0.35 | 0.00 | 0.35 | 0.00 | 0.35 | 0.00 | 0.35 | 0.35 | 0.00 | 0.35 | 0.35 | 0.35 |
| 3 | 0.00 | 0.35 | 0.35 | 0.00 | 0.00 | 0.35 | 0.35 | 0.35 | 0.00 | 0.35 | 0.35 | 0.35 |
| 4 | 0.00 | 0.00 | 0.35 | 0.35 | 0.35 | 0.35 | 0.35 | 0.00 | 0.35 | 0.35 | 0.35 | 0.00 |
| 5 | 0.00 | 0.35 | 0.00 | 0.35 | 0.00 | 0.00 | 0.35 | 0.35 | 0.35 | 0.35 | 0.35 | 0.35 |
| 6 | 0.00 | 0.00 | 0.00 | 0.41 | 0.00 | 0.00 | 0.41 | 0.41 | 0.41 | 0.41 | 0.00 | 0.41 |
| 7 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.00 | 0.33 | 0.00 | 0.33 | 0.33 | 0.00 | 0.33 |
| 8 | 0.00 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.00 | 0.00 |
| 9 | 0.30 | 0.30 | 0.30 | 0.30 | 0.30 | 0.30 | 0.30 | 0.30 | 0.30 | 0.00 | 0.30 | 0.30 |
| 10 | 0.41 | 0.00 | 0.00 | 0.41 | 0.00 | 0.41 | 0.41 | 0.00 | 0.41 | 0.41 | 0.00 | 0.00 |
| 11 | 0.32 | 0.32 | 0.32 | 0.32 | 0.32 | 0.32 | 0.00 | 0.00 | 0.32 | 0.32 | 0.32 | 0.32 |
| 12 | 0.58 | 0.00 | 0.00 | 0.00 | 0.00 | 0.58 | 0.00 | 0.58 | 0.00 | 0.00 | 0.00 | 0.00 |
| 13 | 0.00 | 0.00 | 0.35 | 0.35 | 0.35 | 0.35 | 0.35 | 0.00 | 0.00 | 0.35 | 0.35 | 0.35 |
| 14 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.00 | 0.00 | 0.33 | 0.00 | 0.33 | 0.33 |
| 15 | 0.38 | 0.00 | 0.00 | 0.38 | 0.00 | 0.38 | 0.38 | 0.00 | 0.38 | 0.00 | 0.38 | 0.38 |
# Otras transformaciones posibles con datos de presencia/ausencia
# mc_t <- decostand(mc, 'normalize') #Chord
# mc_t <- decostand(log1p(mc), 'normalize') #Chord
# mc_t <- decostand(mc, 'chi.square') #Chi-square
env <- read_csv(paste0(fuentes_practica, 'matriz-ambiental-', params$estudiante, '.csv'))[, -1]
env %>% estilo_kable(titulo = 'Matriz ambiental', nombres_filas = T, alinear = 'r')
| var1 | var2 | var3 | var4 | |
|---|---|---|---|---|
| 1 | 0.05 | 0.92 | 0.72 | 0.42 |
| 2 | 0.96 | 0.30 | 0.24 | 0.80 |
| 3 | 0.57 | 0.34 | 0.80 | 0.89 |
| 4 | 0.66 | 0.69 | 0.84 | 0.72 |
| 5 | 0.50 | 0.36 | 0.77 | 0.45 |
| 6 | 0.74 | 0.09 | 0.17 | 0.50 |
| 7 | 0.62 | 0.44 | 0.80 | 0.69 |
| 8 | 0.00 | 0.70 | 0.69 | 0.41 |
| 9 | 0.04 | 0.89 | 0.78 | 0.32 |
| 10 | 0.36 | 0.75 | 0.42 | 0.58 |
| 11 | 0.15 | 0.99 | 0.97 | 0.18 |
| 12 | 0.20 | 0.41 | 0.41 | 0.60 |
| 13 | 0.46 | 0.87 | 0.57 | 0.77 |
| 14 | 0.91 | 0.70 | 0.40 | 0.06 |
| 15 | 0.75 | 0.41 | 0.85 | 0.28 |
La matriz ambiental se compone de 4 variables de tipo numérico, conteniendo el valor de cada variable para cada uno de los 15 sitios. La siguiente tabla y el gráfico muestran un resumen de los estadísticos básicos de la matriz ambiental.
estad_basicos <- env %>%
pivot_longer(everything(), names_to = "Variable", values_to = "Valor") %>%
group_by(Variable) %>%
summarise(
Media = mean(Valor),
Mediana = median(Valor),
`Desv. Estándar` = sd(Valor),
Varianza = var(Valor),
`Error Estándar` = sd(Valor) / sqrt(length(Valor)))
estad_basicos %>% estilo_kable(titulo = 'Matriz ambiental', nombres_filas = F, alinear = 'crrrr')
| Variable | Media | Mediana | Desv. Estándar | Varianza | Error Estándar |
|---|---|---|---|---|---|
| var1 | 0.46 | 0.50 | 0.32 | 0.10 | 0.08 |
| var2 | 0.59 | 0.69 | 0.27 | 0.07 | 0.07 |
| var3 | 0.63 | 0.72 | 0.24 | 0.06 | 0.06 |
| var4 | 0.51 | 0.50 | 0.24 | 0.06 | 0.06 |
env %>%
pivot_longer(everything(), names_to = 'Variable', values_to = 'Valor') %>%
group_by(Variable) %>%
ggplot() +
aes(x = Variable, y = Valor, color = Variable, fill = Variable) +
# geom_boxplot(lwd = 0.2) +
geom_violin(alpha = 0.2, width = 0.8, color = "transparent") +
geom_jitter(alpha = 0.6, size = 2, height = 0, width = 0.1) +
geom_boxplot(alpha = 0, width = 0.3, color = "#808080") +
scale_fill_brewer(palette = 'Set1') +
theme_bw() +
theme(legend.position="none")
Las medias calculadas de las variables var1, var2, var3 y var4 son, respectivamente, las siguientes: 0.46, 0.59, 0.63 y 0.51. La variable que con la media más alta fue var3 (0.63), y la más baja la obtuvo la variable var1 (0.46). Por otra parte, la mitad de los sitios midieron menos de 0.5, 0.69, 0.72 y 0.5, para cada una de las variables var1, var2, var3 y var4, respectivamente. Finlamente, la variable con mayor dispersión fue var1 y la de menor dispersión fue var4.
Una verificación importante que debe realizarse es si las matrices de comunidad y ambiental tienen el mismo numero de filas y si las filas se encuentran en el mismo orden (e.g. consistencia entre matrices, donde cada fila en la matriz de comunidad se refiere al mismo sitio en la ambiental, y viceversa). Esto se puede comprobar por medio de los nombres de columnas y, en este caso, tras realizar la correspondiente comprobación, esta condición se cumple, por lo que podemos continuar adelante con los siguientes análisis
A continuación, realizaré análisis de agrupamiento, ordenación y diversidad, basándome en las indicaciones de Borcard, Gillet, y Legendre (2018), reaprovechando el código contenido en Martínez-Batlle (2020).
A continuación, el análisis de agrupamiento propiamente. La parte más importante es generar un árbol, a partir de una matriz de distancias, que haga sentido desde el punto de vista de la comunidad y la distribución de las especies. Primero cargaré paquetes específicos de esta técnica y generaré la matriz de distancias.
mc_d <- vegdist(mc_t, "euc")
A continuación, generaré árboles usando distintos métodos. Explico detalladamente estas técnicas en el repo, y en los vídeos (13 a 16) de la lista mencionada arriba “Ecología Numérica con R” de mi canal.
lista_cl <- list(
cl_single = hclust(mc_d, method = 'single'),
cl_complete = hclust(mc_d, method = 'complete'),
cl_upgma = hclust(mc_d, method = 'average'),
cl_ward = hclust(mc_d, method = 'ward.D2')
)
par(mfrow = c(2,2))
invisible(map(names(lista_cl), function(x) plot(lista_cl[[x]], main = paste0(x, '\n(árbol de evaluación)'), hang = -1)))
par(mfrow = c(1,1))
A continuación, calcularé la distancia y la correlación cofenéticas; esta última, la correlación cofenética,se utiliza como criterio flexible para elegir el método de agrupamiento idóneo, pero no debe usarse de manera estricta. Se supone que el método con la mayor correlación cofenética explica mejor el agrupamiento de la comunidad. Si quieres comprender mejor esta técnica, consulta el vídeo que te referí en el párrafo anterior, así como los libros de referencia. Normalmente, el método UPGMA obtiene la mayor correlación cofenética, pero esto se debe a que su procedimiento de obtención maximiza precisamente dicha métrica. No es recomendable conservar un único método de agrupamiento, normalmente es bueno usar al menos dos. Ward es muchas veces recomendado como método de contraste, por basarse en procedimientos de cálculo muy distintos a los de UPGMA.
map_df(lista_cl, function(x) {
coph_d <- cophenetic(x)
corr <- cor(mc_d, coph_d)
return(corr)
}) %>% t() %>% as.data.frame() %>%
rownames_to_column %>%
mutate(rowname = gsub('cl_', '', rowname)) %>%
setNames(c('Método de agrupamiento', 'Correlación cofenética')) %>%
estilo_kable()
| Método de agrupamiento | Correlación cofenética |
|---|---|
| single | 0.83 |
| complete | 0.81 |
| upgma | 0.86 |
| ward | 0.71 |
Ahora, calcularé las anchuras de silueta, una métrica que ayuda a determinar en cuántos grupos se organiza la comunidad; las anchuras de silueta no deben usarse como método estricto, y sólo debe usarse de forma flexible para informarnos sobre el número máximo de grupos posibles. Considera las siguientes reglas:
# UPGMA
anch_sil_upgma <- calcular_anchuras_siluetas(
mc_orig = mc,
distancias = mc_d,
cluster = lista_cl$cl_upgma)
u_dend_reord <- reorder.hclust(lista_cl$cl_upgma, mc_d)
plot(u_dend_reord, hang = -1, main = 'Método UPGMA\n(árbol de evaluación)')
rect.hclust(
tree = u_dend_reord,
k = anch_sil_upgma$n_grupos_optimo)
resultado_evaluacion_upgma <- evaluar_arbol(u_dend_reord, anch_sil_upgma$n_grupos_optimo)
Tras cortar el árbol, la evaluación practicada concluyó lo siguiente: “Árbol no recomendado para usarse por producir grupos compuestos por dos elementos o menos”
# Ward
anch_sil_ward <- calcular_anchuras_siluetas(
mc_orig = mc,
distancias = mc_d,
cluster = lista_cl$cl_ward)
w_dend_reord <- reorder.hclust(lista_cl$cl_ward, mc_d)
plot(w_dend_reord, hang = -1, main = 'Método Ward\n(árbol de evaluación)')
rect.hclust(
tree = w_dend_reord,
k = anch_sil_ward$n_grupos_optimo)
resultado_evaluacion_ward <- evaluar_arbol(w_dend_reord, anch_sil_ward$n_grupos_optimo)
Tras cortar el árbol, la evaluación practicada concluyó lo siguiente: “Árbol útil para análisis posteriores, siempre que se corte en 2 grupos”.
Una forma alterna de evaluar árboles consiste en usar el remuestreo por bootstrap multiescalar. No me interesa que profundices en ella, sólo presentártela como técnica probabilística para evaluar árboles generados por métodos determinísticos. La técnica es documentada en Borcard, Gillet, y Legendre (2018), de la cual puedes un resumen en este cuaderno y en este vídeo (minuto 51:33). El remuestreo por bootstrap multiescalar valida la robustez de los análisis de agrupamiento tomando múltiples muestras aleatorias de los datos en diferentes tamaños. Este proceso determina qué grupos son consistentemente identificados como clústeres, generando valores de probabilidad aproximadamente insesgados (AU) que son considerados más fiables que las probabilidades de bootstrap tradicionales (BP). Esta técnica ayuda a identificar y confirmar patrones robustos en los datos.
Lo aplicaré primero al árbol generado por el método UPGMA.
# UPGMA
# if(interactive()) dev.new()
cl_pvclust_upgma <-
pvclust(t(mc_t),
method.hclust = "average",
method.dist = "euc",
iseed = 99, # Resultado reproducible
parallel = TRUE, quiet = TRUE)
# Añadir los valores de p
plot(cl_pvclust_upgma, hang = -1, main = 'Método UPGMA bootstrap\n(árbol de evaluación)')
# Añadir rectángulos a los grupos significativos
lines(cl_pvclust_upgma)
pvrect(cl_pvclust_upgma, alpha = 0.90, border = 4)
Lo aplicaré también al árbol generado por el método Ward.
# Ward
# if(interactive()) dev.new()
cl_pvclust_ward <-
pvclust(t(mc_t),
method.hclust = "ward.D2",
method.dist = "euc",
iseed = 99, # Resultado reproducible
parallel = TRUE, quiet = TRUE)
# Añadir los valores de p
plot(cl_pvclust_ward, hang = -1, main = 'Método Ward bootstrap\n(árbol de evaluación)')
# Añadir rectángulos a los grupos significativos
lines(cl_pvclust_ward)
pvrect(cl_pvclust_ward, alpha = 0.91, border = 4)
Basado en lo anterior, elegiré un método de agrupamiento y un número de grupos, y lo exportaré a un archivo que posteriormente podré reaprovechar. La lógica empleada para elegir método de agrupamiento y número de grupos, es la siguiente: si el árbol generado por el método UPGMA no es recomendable (por tener grupos formados 2 o menos elementos), pero Ward sí, se usar el árbol generado por el método Ward y el número de grupos idóneo sugerido por la anchura de silueta. Si UPGMA es recomendable pero Ward no lo es, se usar el árbol generado por el método UPGMA, cortado en el número de grupos sugerido por la anchura de siluetas. Si ambos métodos son recomendables y sugieren el mismo número de grupos, se opta por el arbol generado por el método Ward. Si ambos métodos son recomendables pero sugieren un número diferente de grupos, se elige el método que sugiere menos grupos. Finalmente, si ambos métodos, UPGMA y Ward, resultan ser poco idóneos porque generan grupos muy pequeños (dos o menos elementos), se opta, como último recurso, por elegir el árbol generado por el método Ward cortado en 3 grupos.
grupos_seleccionados <- seleccionar_y_cortar_arbol(
arbol_upgma = lista_cl$cl_upgma, arbol_ward = lista_cl$cl_ward,
resultado_evaluacion_upgma = resultado_evaluacion_upgma,
resultado_evaluacion_ward = resultado_evaluacion_ward)
saveRDS(grupos_seleccionados$resultado,
paste0(fuentes_practica, 'grupos_seleccionados-', params$estudiante,'.RDS'))
El árbol generado por el método UPGMA produce grupos compuestos por dos elementos o menos. Usamos el árbol generado por el método Ward cortado en 2 grupos . El árbol resultante se muestra a continuación:
# Convierte el hclust en dendrograma
dend <- as.dendrogram(grupos_seleccionados$arbol)
# Corta y colorea el dendrograma en k grupos
dend_colored <- color_branches(dend, k=grupos_seleccionados$k)
# Etiqueta los grupos
labels_colors <- labels_colors(dend_colored)
labels(dend_colored) <- paste0(labels(dend_colored), " (",
grupos_seleccionados$resultado[grupos_seleccionados$arbol$order],
")")
# Grafica el dendrograma
# par(mar = c(3, 4, 4, 2) + 0.1) # Ajusta los márgenes
plot(
dend_colored,
main=paste(
'Árbol seleccionado\nMétodo',
grupos_seleccionados$metodo,
'cortado en',
grupos_seleccionados$k, 'grupos'),
xlab = 'Sitios (grupo de pertenencia)')
Apliquemos el análisis de agrupamiento a la matriz ambiental. La clave en este punto es que, si la matriz ambiental presenta patrones parecidos a los de la matriz de comunidad, significa que el agrupamiento utilizado hace sentido entre ambos conjuntos de datos (comunidad y hábitat) de forma consistente. Si ambos conjuntos de datos son consistentes, significa que existe algún grado de asociación, aunque sea sólo una mera asociación estadística.
Agrupar los sitios de muestreo de la matriz ambiental según los grupos previamente definidos.
env_grupos <- env %>%
rownames_to_column('sitios_de_muestreo') %>%
mutate(grupos = as.factor(grupos_seleccionados$resultado)) %>%
pivot_longer(-c(grupos, sitios_de_muestreo), names_to = "variable", values_to = "valor")
Evaluar efectos entre los grupos (“diferencias significativas”). Se utilizan las pruebas estadísticas ANOVA (evalúa homongeneidad de medias) y Kruskal-Wallis (evalúa homogeneidad de medianas). Las tablas están ordenadas en orden ascendente por la columna p_valor_a, que son los p-valores de la prueba ANOVA.
env_grupos_ak <- env_grupos %>%
group_by(variable) %>%
summarise(
p_valor_a = tryCatch(oneway.test(valor ~ grupos)$p.value, error = function(e) NA),
p_valor_k = tryCatch(kruskal.test(valor ~ grupos)$p.value, error = function(e) NA)
) %>%
arrange(p_valor_a)
env_grupos_ak %>% estilo_kable(alinear = 'crr')
| variable | p_valor_a | p_valor_k |
|---|---|---|
| var2 | 0.07 | 0.11 |
| var3 | 0.30 | 0.39 |
| var1 | 0.68 | 0.62 |
| var4 | 0.70 | 0.71 |
Explora tus resultados.
env_grupos %>%
group_by(variable) %>%
ggplot() + aes(x = grupos, y = valor, group = grupos, fill = grupos) +
geom_boxplot(lwd = 0.2) +
scale_fill_brewer(palette = 'Set1') +
theme_bw() +
theme(legend.position="none") +
facet_wrap(~ variable, scales = 'free_y', ncol = 8)
El objetivo de adjuntarle, a la matriz ambiental, el vector de agrupamiento generado a partir de datos de comunidad, consiste en caracterizar ambientalmente los hábitats de los subgrupos diferenciados según su composición. Observa los resultados de las pruebas estadísticas, de los diagramas de caja, y explora tus resultados:
Análisis de preferencia/fidelidad de especies con grupos (clusters), mediante el coeficiente de correlación biserial puntual (phi).
set.seed(9999)
phi <- multipatt(
mc,
grupos_seleccionados$resultado,
func = "r.g",
max.order = 1,
control = how(nperm = 999))
summary(phi)
Multilevel pattern analysis
---------------------------
Association function: r.g
Significance level (alpha): 0.05
Total number of species: 12
Selected number of species: 2
Number of species associated to 1 group: 2
List of species associated to each combination:
Group A #sps. 2
stat p.value
sp03 1.000 0.001 ***
sp05 0.905 0.002 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Tabla de especies que presentaron asociación con grupos por medio de phi, usando umbral de significancia (umbral_alfa).
tabla_phi_sign <- phi$sign
tabla_phi_sign_alfa <- tabla_phi_sign[phi$sign$p.value < umbral_alfa, ]
data.frame(
`Nombre de especie` = rownames(tabla_phi_sign_alfa),
`P-valor` = tabla_phi_sign_alfa$p.value,
`Grupo de asociación` = gsub('s\\.', '', names(tabla_phi_sign_alfa)[tabla_phi_sign_alfa$index]),
check.names = F) %>%
arrange(`Nombre de especie`) %>%
estilo_kable(alinear = 'crr')
| Nombre de especie | P-valor | Grupo de asociación |
|---|---|---|
| sp03 | 0 | A |
| sp05 | 0 | A |
Me basaré en los scripts que comienzan por to_ de este repo, los cuales explico en los vídeos de “Técnicas de ordenación” de la lista de reproducción “Ecología Numérica con R” de mi canal.
pca_mc_t <- rda(mc_t)
summary(pca_mc_t)
Call:
rda(X = mc_t)
Partitioning of variance:
Inertia Proportion
Total 0.3647 1
Unconstrained 0.3647 1
Eigenvalues, and their contribution to the variance
Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11
Eigenvalue 0.08711 0.07232 0.06176 0.04138 0.03694 0.02628 0.01427 0.01225 0.00868 0.002266 0.001216
Proportion Explained 0.23888 0.19833 0.16937 0.11348 0.10130 0.07207 0.03912 0.03360 0.02380 0.006213 0.003334
Cumulative Proportion 0.23888 0.43721 0.60658 0.72006 0.82136 0.89343 0.93255 0.96615 0.98996 0.996170 0.999504
PC12
Eigenvalue 0.0001808
Proportion Explained 0.0004957
Cumulative Proportion 1.0000000
Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions
* General scaling constant of scores: 1.503159
Species scores
PC1 PC2 PC3 PC4 PC5 PC6
sp01 -0.36180 0.0003779 0.19044 0.115317 -0.072784 0.21953
sp02 0.03859 -0.2102061 -0.07232 -0.029765 -0.309445 -0.13125
sp03 0.06625 -0.3708149 -0.07875 -0.085095 0.044868 0.04295
sp04 0.26216 0.0975285 0.26282 -0.003807 -0.034922 -0.04732
sp05 0.06744 -0.3538302 0.03641 -0.106905 0.009449 0.13422
sp06 -0.29189 -0.0468403 0.16414 -0.089150 0.169805 -0.20301
sp07 0.27334 0.1902360 -0.05539 0.002129 0.138968 0.02465
sp08 -0.21938 0.1719524 -0.38232 -0.045473 -0.116248 -0.03812
sp09 0.15193 0.0801882 0.24421 -0.078567 -0.213652 -0.05026
sp10 0.23606 0.0374621 -0.14062 -0.206519 0.053344 0.08132
sp11 0.05740 -0.2355118 -0.04502 0.214424 0.110437 -0.14661
sp12 0.19885 -0.0218691 -0.08346 0.344557 -0.046483 0.03772
Site scores (weighted sums of species scores)
PC1 PC2 PC3 PC4 PC5 PC6
sit1 -0.400636 -0.41514 -0.16809 -0.41080 -0.3157 -0.04846
sit2 -0.034495 -0.18461 -0.61124 0.31495 0.3442 0.91043
sit3 -0.007596 -0.10855 -0.73600 0.13939 0.2264 -0.58329
sit4 0.295998 -0.20006 0.26321 -0.49296 0.5863 -0.21793
sit5 0.411453 0.36118 -0.34645 0.26606 -0.4892 -0.44765
sit6 0.440081 0.78063 -0.26017 0.01407 -0.2341 0.15827
sit7 0.333191 -0.13523 0.16915 -0.06880 -0.4707 0.77600
sit8 0.118074 -0.02619 -0.11408 -0.84165 -0.2190 -0.24415
sit9 -0.107219 -0.19560 0.04735 0.27393 -0.2651 -0.15333
sit10 -0.038826 0.60157 0.61499 -0.42046 0.2292 0.20071
sit11 0.006095 -0.44769 0.29980 0.08632 -0.2433 0.01651
sit12 -1.176808 0.43624 -0.13692 -0.04902 0.1080 0.06170
sit13 0.326735 -0.28058 -0.03951 0.09047 0.8446 -0.02691
sit14 -0.126222 -0.51463 0.44364 0.35977 -0.3427 -0.15652
sit15 -0.039826 0.32864 0.57431 0.73872 0.2413 -0.24537
screeplot(
pca_mc_t,
bstick = TRUE,
npcs = length(pca_mc_t$CA$eig)
)
# Biplot
cleanplot.pca(pca_mc_t, scaling = 1, mar.percent = 0.06, cex.char1 = 0.7)
# Realizar el CA
mc_ca <- cca(mc)
Resumen de análisis de correspondencia.
summary(mc_ca)
Call:
cca(X = mc)
Partitioning of scaled Chi-square:
Inertia Proportion
Total 0.513 1
Unconstrained 0.513 1
Eigenvalues, and their contribution to the scaled Chi-square
Importance of components:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10 CA11
Eigenvalue 0.1189 0.1014 0.09183 0.05761 0.05421 0.03776 0.01955 0.01776 0.01065 0.003096 0.0002778
Proportion Explained 0.2317 0.1977 0.17901 0.11229 0.10567 0.07360 0.03810 0.03462 0.02077 0.006034 0.0005414
Cumulative Proportion 0.2317 0.4294 0.60837 0.72066 0.82634 0.89993 0.93804 0.97266 0.99342 0.999459 1.0000000
Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions
Species scores
CA1 CA2 CA3 CA4 CA5 CA6
sp01 -0.70381 -0.09889 -0.45163 -0.17038 -0.09218 0.362569
sp02 -0.02641 -0.23652 0.34465 0.17323 -0.59775 -0.186595
sp03 0.05511 -0.44171 0.30987 0.08769 0.14818 0.073464
sp04 0.35665 0.11067 -0.39713 0.06155 -0.07478 -0.061477
sp05 0.08002 -0.54373 0.14919 0.16416 0.12118 0.250971
sp06 -0.44158 -0.10414 -0.27342 0.10134 0.27316 -0.357081
sp07 0.32685 0.39612 -0.05208 -0.03883 0.22934 0.041746
sp08 -0.63037 0.67312 0.55267 0.08957 -0.03959 0.003528
sp09 0.16666 0.09719 -0.35003 0.21399 -0.28213 -0.038832
sp10 0.26665 0.21273 0.16049 0.25556 0.19146 0.134196
sp11 0.04742 -0.23688 0.14924 -0.39887 0.12491 -0.256337
sp12 0.24155 0.12385 0.09862 -0.54643 -0.16006 0.086734
Site scores (weighted averages of species scores)
CA1 CA2 CA3 CA4 CA5 CA6
sit1 -1.109027 -0.7437 0.71510 0.9958 -0.31310 -0.04155
sit2 -0.332954 0.1043 1.24732 -1.2098 1.20647 2.30718
sit3 -0.169092 0.4765 1.75596 -0.6005 0.39116 -1.52410
sit4 0.902150 -0.6283 -0.41362 0.9690 1.68626 -0.70636
sit5 0.787733 1.4055 0.68934 -0.4128 -1.40335 -0.91721
sit6 1.020844 2.6519 0.02275 0.1024 -0.41743 0.73232
sit7 0.713555 -0.4167 -0.22753 0.3868 -1.05915 1.95049
sit8 0.143586 0.1794 0.53746 2.1376 -0.06341 -0.41224
sit9 -0.403791 -0.2339 0.07914 -0.4150 -0.58648 -0.19578
sit10 -0.040081 1.0085 -2.47514 1.2245 0.75280 0.35810
sit11 0.035562 -1.1018 -0.28329 -0.1010 -0.64198 0.02016
sit12 -4.980270 1.5451 -0.62571 0.1188 0.86943 0.07961
sit13 0.980912 -0.5954 0.19707 -0.6810 1.96772 -0.29064
sit14 -0.209770 -1.4573 -0.50895 -0.6051 -1.10572 -0.37253
sit15 -0.007538 0.4056 -1.98564 -1.9284 0.04809 -0.84256
Gráfico de sedimentación o screeplot.
# Screeplot
screeplot(mc_ca, bstick = TRUE, npcs = length(mc_ca$CA$eig))
Representación del biplot.
# Biplot
plot(mc_ca,
scaling = 1,
main = "Análisis de correspondencia, escalamiento 1"
)
A continuación, el análisis de ordenación propiamente. La parte más importante es el entrenamiento: la función train del paquete caret, contenida en la función my_train, simplifica la selección de variables. Lo más importante: prueba con todas las variables primero, observa las variables que recomienda el modelo final (print_my_train(mod)) y ensaya varias combinaciones de subconjuntos de variables.
mc_t_ren <- mc_t %>%
rename_all(~ paste('ESPECIE', .x))
env_spp <- env %>% bind_cols(mc_t_ren)
spp <- paste0('`', grep('^ESPECIE', colnames(env_spp), value = T), '`', collapse = ' + ')
my_formula <- as.formula(paste(spp, '~ .'))
set.seed(1); mod <- my_train(
formula = my_formula,
# preproceso = 'scale',
data = env_spp,
num_variables = 3:4)
print_my_train(mod)
$resumen_variables
Subset selection object
4 Variables (and intercept)
Forced in Forced out
var1 FALSE FALSE
var2 FALSE FALSE
var3 FALSE FALSE
var4 FALSE FALSE
1 subsets of each size up to 3
Selection Algorithm: 'sequential replacement'
var1 var2 var3 var4
1 ( 1 ) " " " " "*" " "
2 ( 1 ) " " "*" "*" " "
3 ( 1 ) "*" "*" "*" " "
$resultados_nvmax
nvmax RMSE Rsquared MAE RMSESD RsquaredSD MAESD
1 3 0.3621768 0.3274227 0.2737789 0.1078299 0.4043360 0.05911840
2 4 0.3912243 0.2715631 0.3067771 0.1095387 0.3973059 0.04482865
$mejor_ajuste
nvmax
1 3
(covar <- grep(
pattern = '\\(Intercept\\)',
x = names(coef(mod$finalModel,unlist(mod$bestTune))),
invert = T, value = T))
[1] "var1" "var2" "var3"
rda_mc_t <- rda(mc_t_ren %>% rename_all(~ gsub('^ESPECIE ', '', .)) ~ .,
env %>% select_at(all_of(gsub('\\`', '', covar))), scale = T)
A continuación, el resumen del análisis de redundancia.
summary(rda_mc_t)
Call:
rda(formula = mc_t_ren %>% rename_all(~gsub("^ESPECIE ", "", .)) ~ var1 + var2 + var3, data = env %>% select_at(all_of(gsub("\\`", "", covar))), scale = T)
Partitioning of correlations:
Inertia Proportion
Total 12.000 1.0000
Constrained 3.302 0.2752
Unconstrained 8.698 0.7248
Eigenvalues, and their contribution to the correlations
Importance of components:
RDA1 RDA2 RDA3 PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
Eigenvalue 1.9631 1.00273 0.33649 2.4945 2.1649 1.2662 0.92060 0.59610 0.51854 0.3924 0.24687 0.081650
Proportion Explained 0.1636 0.08356 0.02804 0.2079 0.1804 0.1055 0.07672 0.04968 0.04321 0.0327 0.02057 0.006804
Cumulative Proportion 0.1636 0.24715 0.27519 0.4831 0.6635 0.7690 0.84570 0.89538 0.93859 0.9713 0.99186 0.998667
PC10 PC11
Eigenvalue 0.0112046 0.0047924
Proportion Explained 0.0009337 0.0003994
Cumulative Proportion 0.9996006 1.0000000
Accumulated constrained eigenvalues
Importance of components:
RDA1 RDA2 RDA3
Eigenvalue 1.9631 1.0027 0.3365
Proportion Explained 0.5945 0.3036 0.1019
Cumulative Proportion 0.5945 0.8981 1.0000
Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions
* General scaling constant of scores: 3.600206
Species scores
RDA1 RDA2 RDA3 PC1 PC2 PC3
sp01 0.02233 -0.19725 -0.214654 -0.49209 -0.6099 -0.27616
sp02 0.19407 -0.19146 0.484987 -0.11312 0.3868 -0.64242
sp03 0.57826 -0.07897 0.001634 -0.27793 0.7126 -0.01761
sp04 0.28678 0.10262 -0.013890 0.87577 -0.2854 -0.15229
sp05 0.63310 -0.21485 -0.178746 -0.08989 0.5323 -0.17307
sp06 0.24287 -0.49302 0.044248 -0.38293 -0.5291 0.43226
sp07 -0.22370 0.41694 0.082927 0.62464 0.1692 0.42297
sp08 -0.84854 -0.23766 0.071236 -0.38047 0.2243 -0.03025
sp09 0.18266 -0.03336 0.154699 0.70127 -0.3373 -0.46724
sp10 -0.10225 0.03597 0.004855 0.52308 0.6208 0.32306
sp11 0.62497 0.30823 0.108810 -0.37377 0.2037 0.10604
sp12 0.03647 0.61137 -0.017530 0.02135 0.2112 -0.35296
Site scores (weighted sums of species scores)
RDA1 RDA2 RDA3 PC1 PC2 PC3
row1 0.41522 -2.1133 0.827471 -0.80556 0.556417 -0.40180
row2 -0.28045 1.0844 -2.527745 -1.05660 1.295665 0.28504
row3 -0.50612 0.5630 2.528474 -1.11851 0.984571 0.97367
row4 1.06307 -0.1294 -0.635410 0.49831 0.009481 1.52098
row5 -0.83742 1.7519 3.120194 0.71548 0.265308 -0.63986
row6 -1.72332 1.7978 -0.004235 1.47362 0.216896 -0.52610
row7 0.36971 0.5879 -0.023833 0.50695 0.378282 -1.19562
row8 -0.12314 -1.4037 1.608370 0.83436 0.823529 0.02357
row9 0.44289 -0.1350 0.807652 -0.12164 -0.008468 -0.79289
row10 -0.68151 -0.5127 -1.221340 1.44075 -1.341318 0.88662
row11 1.24890 -0.4576 0.223684 -0.07623 -0.005199 -0.70086
row12 -1.87077 -2.9071 -2.621653 -1.84332 -1.477299 0.48697
row13 0.91163 1.0993 -1.657859 0.30776 0.827001 1.49917
row14 1.47582 -0.5338 0.308627 -0.68262 -0.484011 -1.69531
row15 0.09549 1.3082 -0.732397 -0.07274 -2.040856 0.27643
Site constraints (linear combinations of constraining variables)
RDA1 RDA2 RDA3 PC1 PC2 PC3
row1 0.34115 -1.3282 0.14057 -0.80556 0.556417 -0.40180
row2 -0.41630 1.0095 -1.47009 -1.05660 1.295665 0.28504
row3 -0.45663 0.8635 1.12532 -1.11851 0.984571 0.97367
row4 1.15501 0.7586 0.34593 0.49831 0.009481 1.52098
row5 -0.57968 0.6010 1.04918 0.71548 0.265308 -0.63986
row6 -1.89621 0.5432 -0.93284 1.47362 0.216896 -0.52610
row7 0.04192 0.8772 0.82077 0.50695 0.378282 -1.19562
row8 -0.62355 -1.2364 0.59434 0.83436 0.823529 0.02357
row9 0.32716 -1.2163 0.40961 -0.12164 -0.008468 -0.79289
row10 0.04517 -0.8259 -1.06007 1.44075 -1.341318 0.88662
row11 1.22353 -0.7402 0.77004 -0.07623 -0.005199 -0.70086
row12 -1.66263 -0.8715 -0.08379 -1.84332 -1.477299 0.48697
row13 0.95189 -0.4368 -0.87536 0.30776 0.827001 1.49917
row14 1.18111 0.6558 -1.71901 -0.68262 -0.484011 -1.69531
row15 0.36804 1.3465 0.88539 -0.07274 -2.040856 0.27643
Biplot scores for constraining variables
RDA1 RDA2 RDA3 PC1 PC2 PC3
var1 0.02854 0.91152 -0.4103 0 0 0
var2 0.73492 -0.67797 0.0160 0 0 0
var3 0.51923 -0.06324 0.8523 0 0 0
La varianza ajustada explicada por el modelo.
RsquareAdj(rda_mc_t)$adj.r.squared
[1] 0.07751792
Y el factor de inflación de la varianza.
vif.cca(rda_mc_t)
var1 var2 var3
1.624959 1.701148 1.278334
Represento el gráfico triplot.
# Triplot
escalado <- 1
plot(rda_mc_t,
scaling = escalado,
display = c("sp", "lc", "cn"),
main = paste("Triplot de RDA especies ~ variables, escalamiento", escalado)
)
rda_mc_t_sc1 <- scores(rda_mc_t,
choices = 1:2,
scaling = escalado,
display = "sp"
)
# text(mi_fam_t_rda, "species", col="red", cex=0.8, scaling=escalado)
arrows(0, 0,
rda_mc_t_sc1[, 1] * 0.9,
rda_mc_t_sc1[, 2] * 0.9,
length = 0,
lty = 1,
col = "red"
)
Me basaré en los scripts que comienzan por di_ de este repo, los cuales explico en los vídeos de “Análisis de diversidad” (vídeos 19 y 20) de la lista de reproducción “Ecología Numérica con R” de mi canal. Dichos vídeos tienen aplicaciones ligeramente diferentes, pues los datos fuente usados en ellos son de abundancia, mientras que los tuyos son de presencia/ausencia.
La principal desventaja de trabajar con registros de presencia, es que la mayoría de los índices de diversidad alpha fueron diseñados originalmente para calcularse a partir de datos de abundancia. Sin embargo, la riqueza de especies, que es el número \(q=0\) de Hill (\(=N_0\) en las columnas que produce la función alpha_div) es un buen proxy sobre la diversidad, y nos ayudará a comparar sitios.
Además de la columna N0 del objeto que generaré en el bloque siguiente, verás que la función alpha_div genera otras columnas; son índices pensados para datos de abundancia, que en este caso no usaremos, pero los muestro para que tengas una visión completa del análisis de diversidad con índices que podría serte de utilidad en el futuro.
Por otra parte, afortunadamente, los métodos de estimación de riqueza de Chao, y los de diversidad beta (al final de esta sección), aprovechan sustancialmente los registros de presencia/ausencia para realizar estimaciones consistentes y fiables.
Una nota adicional. En el análisis de diversidad, es útil (no imprescindible) disponer de un análisis clúster (agrupamiento) básico. Este te servirá para comparar la riqueza observada y la esperada entre hábitats. Por esta razón, combinamos análisis de diversidad con agrupamiento. Sin embargo, si el análisis de agrupamiento generó grupos de dos o menos elementos, dicha comparación no será realizable.
indices <- alpha_div(mc) %>%
mutate(sitio = rownames(.)) %>%
relocate(sitio, .before = everything())
El objeto mc es la matriz de comunidad de presecia/ausencia. La función alpha_div es un “envoltorio” generado por mí para calcular múltiples índices de diversidad y estimaciones, basada en las funciones de los paquetes SpadeR y iNEXT. Si usásemos datos de abundancia, los índices que calcula la función “alpha_div” serían útiles, pero con registros de presencia/ausencia, como es nuestro caso, sólo la columna N0 (riqueza) nos aportará algún resultado con sentido.
indices %>%
kable(booktabs=T) %>%
kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
gsub(' NA |NaN ', '', .) #Lista de especies
| sitio | N0 | H | Hb2 | N1 | N1b2 | N2 | J | E10 | E20 |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 2 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 3 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 4 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 5 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 6 | 6 | 1.791759 | 2.584963 | 6 | 6 | 6 | 1 | 1 | 1 |
| 7 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 8 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 9 | 11 | 2.397895 | 3.459432 | 11 | 11 | 11 | 1 | 1 | 1 |
| 10 | 6 | 1.791759 | 2.584963 | 6 | 6 | 6 | 1 | 1 | 1 |
| 11 | 10 | 2.302585 | 3.321928 | 10 | 10 | 10 | 1 | 1 | 1 |
| 12 | 3 | 1.098612 | 1.584963 | 3 | 3 | 3 | 1 | 1 | 1 |
| 13 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 14 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 15 | 7 | 1.945910 | 2.807355 | 7 | 7 | 7 | 1 | 1 | 1 |
Los sitios ordenados en función de su riqueza:
indices %>%
arrange(desc(N0)) %>%
kable(booktabs=T) %>%
kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
gsub(' NA |NaN ', '', .) #Lista de especies
| sitio | N0 | H | Hb2 | N1 | N1b2 | N2 | J | E10 | E20 |
|---|---|---|---|---|---|---|---|---|---|
| 9 | 11 | 2.397895 | 3.459432 | 11 | 11 | 11 | 1 | 1 | 1 |
| 11 | 10 | 2.302585 | 3.321928 | 10 | 10 | 10 | 1 | 1 | 1 |
| 1 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 7 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 8 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 14 | 9 | 2.197225 | 3.169925 | 9 | 9 | 9 | 1 | 1 | 1 |
| 2 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 3 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 4 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 5 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 13 | 8 | 2.079442 | 3.000000 | 8 | 8 | 8 | 1 | 1 | 1 |
| 15 | 7 | 1.945910 | 2.807355 | 7 | 7 | 7 | 1 | 1 | 1 |
| 6 | 6 | 1.791759 | 2.584963 | 6 | 6 | 6 | 1 | 1 | 1 |
| 10 | 6 | 1.791759 | 2.584963 | 6 | 6 | 6 | 1 | 1 | 1 |
| 12 | 3 | 1.098612 | 1.584963 | 3 | 3 | 3 | 1 | 1 | 1 |
En el bloque siguiente, represento gráficamente la correlación entre la riqueza y las variables ambientales mediante un panel de gráficos, que suele llamarse también “matriz de correlación”, expresada gráficamente. Si usases índices de diversidad, como el de Shannon o los números de Hill, también deberías incluirlos en el gráfico; nota que en este ejemplo, sólo uso la riqueza (la función select(N0) se encarga de conservar sólo la riqueza). Esto es lo que debes saber sobre el panel:
Presta atención a la primera columna y la primera fila de la matriz, que muestra cómo se correlaciona N0 con las variables ambientales que elijas.
La diagonal contiene gráficos de línea que muestra la densidad de la variable en cuestión.
Los gráficos del “triángulo superior”, y que contienen el patrón Corr: ####, muestran el valor del coeficiente de correlación de Pearson (\(r\)) entre las variables intersectadas. Si existe un \(|r|\) elevado (es decir, si es muy cercano a -1 o a 1) y la prueba de producto-momento es significativa (si hay uno o varios asteriscos, o un punto, lo es), entonces toma nota de que dicha variable se asocia estadísticamente con la riqueza. Si \(r\) es negativo, la relación es inversa (cuando aumenta la variable, disminuye la riqueza, y viceversa); si es positivo, la relación es directa (cuando aumenta la variable, aumenta también la riqueza).
En el “triángulo inferior”, que es un espejo del superior, se sitúan los gráficos de dispersión de las variables intersectadas. Si los puntos siguen un patrón de distribución formando una elipse imaginaria (organizados en torno a una línea recta imaginaria inclinada), entonces existe correlación.
bind_cols(indices %>% select(N0), env %>%
rename_with(.fn = ~ paste0('AMB_', .))) %>%
ggpairs(
labeller = label_wrap_gen(width=10),
upper = list(continuous = wrap("cor", size = 3))) +
theme(text = element_text(size = 10))
“Completitud”, en porcentajes, según distintos estimadores. Con un 80% de completitud, se considera en general una muestra representativa. Sin embargo, este umbral de 80% no debe tomarse de forma estricta. Sobre todo porque existen métodos refinados que mejoran las estimaciones
riqueza_estimaciones <- data.frame(specpool(mc) %>% select(-matches('.se$'))) %>%
select(`Riqueza observada` = Species,
`Número de sitios` = n,
`Estimación por Chao (clásico)` = chao,
`Estimación por jackknife de primer orden` = jack1,
`Estimación por jackknife de segundo orden` = jack2,
`Estimación por bootstrap` = boot) %>%
pivot_longer(cols = everything(), names_to = 'Variable', values_to = 'Valor') %>%
mutate(`Cobertura (%)` = Valor / (filter(., Variable == "Riqueza observada") %>% pull(Valor)) * 100) %>%
mutate(`Cobertura (%)` = ifelse(Variable %in% c('Riqueza observada', 'Número de sitios'), NA, `Cobertura (%)`))
riqueza_estimaciones %>% estilo_kable(alinear = 'lrr')
| Variable | Valor | Cobertura (%) |
|---|---|---|
| Riqueza observada | 12 | |
| Número de sitios | 15 | |
| Estimación por Chao (clásico) | 12 | 100 |
| Estimación por jackknife de primer orden | 12 | 100 |
| Estimación por jackknife de segundo orden | 12 | 100 |
| Estimación por bootstrap | 12 | 100 |
# Bug no resuelto:
# Error in if (var_mle > 0) { : valor ausente donde TRUE/FALSE es necesario
# Varios intentos frustrados por lograr que funcione. Entiendo que el problema
# está en el número de doubletons (la matriz no tiene), pero no logré mejorar
# la función interna SpecInciHomo para solucionarlo. La versión de SpadeR usada
# en la aplicación Shiny https://chao.shinyapps.io/SpadeR/, no es la misma que
# la que se encuentra en GitHub ni en el CRAN, pues esa no tiene bug.
df_spader <- data.frame(V1 = as.integer(c(nrow(mc), colSums(mc))))
# También se puede crear con esta línea:
# df_spader <- structure(
# list(V1 = c(15, 8, 9, 10, 9, 8, 9, 9, 8, 8, 6, 12, 9)),
# class = "data.frame", row.names = c(NA, -13L))
df_spader
# V1
# 15
# 8
# 9
# 10
# 9
# 8
# 9
# 9
# 8
# 8
# 6
# 12
# 9
ChaoSpecies(df_spader, datatype = 'incidence_freq',
k = min(df_spader$V1), conf=0.95)
# Error in if (var_mle > 0) { : valor ausente donde TRUE/FALSE es necesario
# ENG: Error in if (var_mle > 0) { : missing value where TRUE/FALSE needed
Graficaré la curva de acumulación de especies.
mc_general <- mc %>%
summarise_all(sum) %>%
mutate(N = nrow(mc)) %>%
relocate(N, .before = 1) %>%
data.frame
nasin_raref <- iNEXT::iNEXT(
x = t(mc_general),
q=0,
knots = 2000,
datatype = 'incidence_freq')
acumulacion_especies <- iNEXT::ggiNEXT(nasin_raref, type=1) +
theme_bw() +
theme(
text = element_text(size = 20),
panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
) +
ylab('Riqueza de especies') +
xlab('Número de sitios') +
scale_y_continuous(breaks = seq(0, 80, length.out = 9)) +
scale_color_manual(values = brewer.pal(8, 'Set2')) +
scale_fill_manual(values = brewer.pal(8, 'Set2'))
acumulacion_especies
Ahora según los grupos previamente seleccionados en el análisis de agrupamiento.
grupos_seleccionados <- readRDS(paste0(
fuentes_practica, 'grupos_seleccionados-',
params$estudiante, '.RDS'))
mc_grupos <- mc %>%
mutate(g = grupos_seleccionados) %>%
group_by(g) %>%
summarise_all(sum) %>%
select(-g) %>%
mutate(N = nrow(mc)) %>%
relocate(N, .before = 1) %>%
data.frame
nasin_raref_general <- iNEXT::iNEXT(
x = t(mc_grupos),
q=0,
knots = 400,
datatype = 'incidence_freq')
acumulacion_especies_grupos <- iNEXT::ggiNEXT(nasin_raref_general, type=1) +
theme_bw() +
theme(
text = element_text(size = 20),
panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
) +
ylab('Riqueza de especies') +
xlab('Número de sitios') +
scale_y_continuous(breaks = seq(0, 80, length.out = 9)) +
scale_color_manual(values = brewer.pal(8, 'Set2')) +
scale_fill_manual(values = brewer.pal(8, 'Set2'))
acumulacion_especies_grupos
determinar_contrib_local_y_especie(
mc = mc,
alpha = 0.05,
nperm = 9999,
metodo = 'sorensen')
## $betadiv
## $beta
## SStotal BDtotal
## 2.6200609 0.1871472
##
## $SCBD
## [1] NA
##
## $LCBD
## [1] 0.05822156 0.06974966 0.06656907 0.05190575 0.06325906 0.08934245 0.05035117 0.04897587 0.02940812 0.08071253
## [11] 0.03715550 0.17737369 0.05516652 0.05387180 0.06793727
##
## $p.LCBD
## [1] 0.5766 0.3987 0.4416 0.6838 0.4941 0.1749 0.7030 0.7274 0.9619 0.2803 0.8970 0.0009 0.6351 0.6569 0.4271
##
## $p.adj
## [1] 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.0135 1.0000 1.0000 1.0000
##
## $method
## [1] "sorensen" "sqrt.D=FALSE"
##
## $note
## [1] "Info -- D is Euclidean because beta.div outputs D[jk] = sqrt(1-S[jk])"
## [2] "For this D functions, use beta.div with option sqrt.D=FALSE"
##
## $D
## [1] NA
##
## attr(,"class")
## [1] "beta.div"
##
## $especies_contribuyen_betadiv
## [1] NA
##
## $sitios_contribuyen_betadiv
## [1] "12"
##
## $valor_de_ajustado_lcbd
## [1] 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.0135 1.0000 1.0000 1.0000
##
## $sitios_contribuyen_betadiv_ajustado
## [1] "12"